Discrete stochastic modeling of calcium channel dynamics 
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We propose a simple discrete stochastic model for calcium dynamics in living cells. Specifically, 
the calcium concentration distribution is assumed to give rise to a set of probabilities for the open- 
ing/closing of channels which release calcium thereby changing those probabilities. We study this 
model in one dimension, analytically in the mean-field limit of large number of channels per site N, 
and numerically for small N. As the number of channels per site is increased, the transition from a 
non-propagating region of activity to a propagating one changes in nature from one described by di- 
rected percolation to that of deterministic depinning in a spatially discrete system. Also, for a small 
number of channels a propagating calcium wave can leave behind a novel fluctuation-driven state, 
in a parameter range where the limiting deterministic model exhibits only single pulse propagation. 

PACS: 87.16. Xa,05.40.-a,82.20.Mj 

It has become clear that the intracellular nonlinear dy- 
namics of calcium plays a crucial role in many biological 
processes |l|. The nonlinearity of this problem is due 
to the fact that there exist calcium stores inside the cell 
which can be released via the opening of channels which 
themselves have calcium-dependent kinetics. Typically, 
these processes are modeled using a set of coupled equa- 
tions for the calcium concentration (the diffusion equa- 
tion with sources and sinks) and for the relevant chan- 
nels; the latter is often described by a rate equation for 
the fraction of open channels per unit of area. More 
elaborate models take into account the discrete nature of 
these channels, their spatial clustering, and fluctuations 
in the process of their opening and closing [g|| . 

In this paper, we will propose and analyze a set of mod- 
els which operate just with the channel dynamics alone. 
The justification for this is that the calcium field equi- 
librates quickly, with a diffusion time of perhaps 0.1s, 
as compared to the channel transition times, perhaps on 
the order of Is for activation of a subunit to several sec- 
onds for its deactivation. One can then imagine solving 
for the quasi-stationary calcium concentration and there- 
after using it to determine the conditional probabilities 
of channel opening or closing. In a subsequent paper 
H , we will show how this can be done in detail starting 
from a specific fully-coupled model (the DeYoung-Keizer- 
model [|]||); here, we will make reasonable assumptions 
for these probabilities and study the resulting stochastic 
model in a one dimensional geometry. 

For specificity, we will focus on systems that have IP3 
(inositol 1,4,5- trisphosphate) channels. Each of these 
channels consists of a number of subunits. Here we as- 
sume that h subunits have to be activated for the channel 
to be open; experiments indicate that ft = 3fl, A sub- 
unit is activated when IP3 ion is bound to its correspond- 
ing domain and Ca 2+ is bound to its activating domain 
and not bound to its inhibiting site. The characteristic 
time of binding and unbinding of IP3 is typically so fast 
(more than 20 times faster than other binding steps ||), 



that we can assume local balance of active/passive chan- 
nels maintained at all times. Furthermore, we assume 
that the channels are spatially organized into clusters 
HII , with a fixed number of channels N per cluster and 
a fixed inter-cluster distance. 

Our model is as follows. We introduce two stochas- 
tic variables for each channel cluster: rii, the number 
of activated subunits, and mj, the number of inhibited 
subunits. At every time step, the number of activated 
subunits rii at site i is changed due to three stochastic 
processes; activation of additional subunits by binding 
available Ca 2+ to their activation domains, de-activation 
by unbinding Ca 2+ from active subunits, and inhibition 
by binding available Ca 2+ to their inhibition domains. 
We take these transition rates to depend on the num- 
ber of open channels at site i, Ci, and on the number 
of open channels at the nearest neighboring sites i ± 1, 
Ci±\. Similarly, there will be binding and unbinding to 
the inhibitory domain, changing m.;. We denote by P 0(11 
the probability to activate/inhibit a subunit per number 
of open channels at the same site (0) or the neighbor- 
ing site (1). To compute the actual probabilities, we 
need to multiply these by the number of open channels. 
Here, we use the simple expedient of taking this to equal 
n'l/hN!}~ 1 where the total number of subunits N s = hN; 
this is easily shown to be the expected number of open 
channels for large enough N. This approach allows us 
to avoid keeping explicit account of each of the indepen- 
dent subunits. Also, we let p d be the deactivation and 
deinhibition probabilities which are c independent. 

Let us define the total probabilities p ± = p + 2p 1 
and the "diffusion constant" a — Pi /(pq + 2p 1 ). We 
also denote Ci(t) = (1 — 2a)a(t) + aci_i(t) + ac i+1 (t), 
which mimics the amount of calcium at site i due to open 
channels at sites i, i±l. Our model explicitly consists of 
the following coupled stochastic processes, rii is updated 

n l (t + At) = n l (t) + A+-A--S+ (1) 

where A+ is a random integer number drawn from the 



binomial distribution B(A+, N s — rii(t) — rrii(t) , p + Ci{t)) , 
A^ is drawn from B(A~, rii(t) , p~ Ci(t)) , and 5£ is drawn 
from B(5n ,rii(t),pj). The equation for rrii reads 



rrii(t + At) = mi {t) + A r H 



(2) 



where A+ is drawn from B(A^ l ,N s — irii(t) , p~~ Ci(t)) , 
and 5+ is drawn from B(5^, mi(t),pj). We do not allow 
for transitions from the inhibited state to the activated 
state. In all these formulas, B(x, y,p) = {^)p x (l — p) y ~ x . 
Note that the probability that IP 3 is bound is included 
by rescaling the number of subunits. 

As a first step, we consider a simplified version of the 
channel dynamics with the inhibition process excluded 
(all p~—0), i.e. a subunit is activated whenever Ca 2+ is 
attached to its activating site. Thus we take rrii = 0, 
and arrive at the one-variable model for the number of 
activated subunits rij. Let us first focus on fairly small 
N s . Examples of the stochastic dynamics for several val- 
ues of parameters are shown in Figure [D. At small a, an 
initial seed almost always ultimately dies giving rise to 
so-called abortive calcium waves. At larger values of a 
the region of activated channels typically expands at a 
finite rate. This transition mirrors what has been seen in 
many experimental systems 0. 
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FIG. 1. Evolution of an initial seed of 5 clusters of open 
channels in the middle of the lattice for p + N s /h — 1, 
p+ = 0.2, h = 3, and N 3 = 3, a = 0.1 (a), N s = 3, a = 0.25 
(b), and N s — 30, a = 0.25 (c). The horizontal axis is spa- 
tial coordinate (100 sites), and the vertical axis is time (1000 
iterations) . 

As is well known for statistical models such as the con- 
tact process Eu|, the critical value of a can be accu- 
rately determined by computing the distribution of sur- 
vival times II(i) for the activation process started from a 
single active site. For a < a c , the distribution falls expo- 
nentially at large t as the wave of activation eventually 
dies out. On the contrary, at a > a c , LT(t) asymptotically 
reaches a constant value 11^ , since a non-zero fraction of 
runs produce ever-expanding active regions. At a — a c , 
the distribution function exhibits a power-law asymptotic 
behavior with the slope determined by the universality 
class of the underlying stochastic process. Our data (not 
shown) indicate that a c is inversely proportional to the 
number of subunits per site N s . We have checked that 
our data is in the directed percolation (DP) |11[| class. For 
example, in Fig. gwe show II(£) of a cluster of open chan- 
nels at the critical value of a c for h — 3, N s — 10 and 



7 = 0.1. The power-law dependence is consistent with 
DP prediction of LT(t) oc t~ - 159 . This is perhaps not too 
surprising. According to the Janssen-Grassberger DP 
conjecture p2|, any spatio-temporal stochastic process 
with short range interactions, fluctuating active phase 
and unique non-fluctuating (absorbing) state, single or- 
der parameter and no additional symmetries, should be- 
long to the DP class. This result does open up the excit- 
ing possibility that intracellular calcium dynamics could 
be an experimental realization of the DP process. 




FIG. 2. The distribution of survival time for the stochastic 
model with h — 3, 7 = 0.1, N s = 10, and a — a c — 0.359. 
Dashed line indicates the power-law scaling oc i -0159 . 

Figure [j](c) shows the opposite limit where the dynam- 
ics becomes almost deterministic. If we take iV s — > 00 
and fix pN s /h — > P, we can use a mean-field description 
in terms of the fraction of activated subunits pi = n,i/N s , 



((1 - 2a) p^ + apt, + 



cvtiX 1 



Pi) - jp l . (3) 



and where we rescaled time t' = Pt/ At and introduced 
7 = Pd/P- For all h > 2, if 7 < j cr Eq.(||) the system 
possesses two stable uniform solutions, p = and p = pa 
and one unstable solution p u , where po,u are real roots 
of the algebraic equation /9 /l_1 (l — p) — 7. The front is 
a solution connecting these two stable fixed points; it is 
easy to show that this front has a unique propagation 
velocity. 

For small a, the discreteness of our spatial lattice 
causes the front to become pinned, as the probability of 
activating subunits at the neighboring site O(ciPq) be- 
comes smaller than the threshold value for excitation 
probability 0(p u ). The stationary front solution is de- 
scribed by the recurrence relation, 



(1 - 2a) p 1 ? + 



api-i 



aPi+i 



IPi 
I -Pi 



(4) 



The bifurcation line which separates pinned and moving 
fronts, can be found in the limit of small a by using the 
ideas of ref. [pL3|| . Indeed, in this limit, the values of pi 
quickly (as a 1 ) approach and po away from the front 
at i — ► ±00, respectively. We can thus replace pi by 
po and everywhere to the left and to the right of the 
front position except for p± at the two sites nearest to 



the front, i — 1 and i + 1. Solving the resulting set of two 
algebraic equations up to a 2 , one can obtain the values 
of p±. At any 7, there is a critical value of a m at which 
the real solution p± vanishes. The family of these values 
a m forms the bifurcation line for front pinning in (7, a) 
plane. At large a, discreteness of the mean field model 
(||) becomes insignificant, and (||) can be replaced by its 
continuum limit 



d tP =(p h -adlp h )(l-p)- 1P . 



(5) 



which of course has no front pinning. Instead, a can 
be scaled out and there is specific value of 7 at which 
the system goes from forward to backward propagating 
fronts. Figure H shows the phase diagram of the mean 
field equation (||) for h = 3. All the data (except pos- 
sibly at the non-generic case 7 = 0) are consistent with 
expected Q (a — am) 1 / 2 scaling. 
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FIG. 3. Phase diagram of the mean field equation (pi) for 
h — 3: a bifurcation line separating forward-propagating 
fronts from pinning region; b same for backward propagat- 
ing fronts; c small-a approximation of pinning line; d line 
7 = 4/27 separating the region of non-existence of the ex- 
cited state; e Maxwell line 7 = 0.138... separating forward 
and backward front propagation in the continuum limit 

How does one get from DP behavior to deterministic 
pinning/dcpinning? To investigate this issue, we have 
performed simulations for the front speed as a function 
a at various finite values of N s ; with the results given in 
Fig. |j. At large N s , the velocity approaches the mean 
field prediction as long as a > a m . Close to critical value 
a m , the velocity deviates from the mean-field dependence 
V oc (a — am) 1 ' 2 because of thermally activated "creep"; 
fluctuations let the front to overcome potential barriers 
associated with finite site separation, and lead to expo- 
nentially slow front propagation (see, e.g., p4|). Directed 
percolation regime is not observed at large N s since the 
DP critical value a c is less than a m . At smaller N s , the 
relative magnitude of the fluctuations grows, and the DP 
threshold value a c exceeds a m . Now, the front pinning 
is determined by fluctuations rather than discreteness, 
and the critical state exhibits the properties of directed 
percolation. 



Now we return to the full two- variable stochastic model 
which describes both activation and inhibition. Since the 
probability of Ca 2+ binding to the inhibition domain is 
typically much smaller than those for the activation do- 
main, the inhibitor dynamics is slow. In the mean-field 
limit N s — ► 00, this model is similar to the FitzHugh- 
Nagumo model often used to describe waves propagat- 
ing in excitable systems. One therefore expects that for 
a certain range of binding/unbinding probabilities, the 
model gives rise to pulse propagation; that is, once the 
wave passes, the system goes into a state dominated by 
inhibition from which it slowly recovers as the inhibitory 
domains slowly unbind. This is indeed what we find for 
large enough N Sl as shown in Fig. |§(a). Behind the 
pair of outgoing pulses, the channels stay refractory for 
a certain time 0(l/p~) and then return to the quiescent 
state. 
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FIG. 4. The average front speed as a function of a for 
stochastic model at h — 3, 7 = 0.1, p + N s /h = 1, pj = 0.1, 
and different values of N„ . Solid line indicates the mean-field 
limit N s — » 00. 

However, we find that having only a modest number 
of channels N leads to fluctuations which strongly affect 
the spatio-temporal behavior of the model. In fact, a new 
dynamical state is formed behind the outgoing fronts, a 
state which remains active at all subsequent times (see 
Fig.||,b). This state is catalyzed by backfiring, i.e. the 
creation of oppositely propagating waves behind a mov- 
ing front. In the deterministic limit of our model, this 
cannot occur as the system is completely refractory once 
the front has passed. At finite N however, propagation of 
the front does not lead to the activation and subsequent 
inhibition of all the channels. Instead, a finite number 
of these remain inactivated, providing a supply of active 
elements that can still support wave propagation. There 
exist more complicated deterministic models |15| ] , such as 
one proposed for CO oxidation on single crystal surfaces 
p6| , which also appear to have pulse-induced backfiring. 
There, however, this effect is due to the loss of pulse sta- 
bility which occurs due to the rather complex non-linear 
dynamics of the inhibitory field. Here, it is the fluctua- 
tions which allow for this phenomenon. 

We have checked that this backfiring-induced state oc- 



curs as well in more realistic and more complex models 
which solve for the calcium concentration together with 
the channel dynamics. Again, the mechanism appears 
to be the lack of complete inhibition in the wake of the 
propagating pulse. Hence, our result that one should find 
this behavior in intracellular calcium dynamics is not an 
artifact of any of the simplifying assumptions used here. 
Also, this state persists when the model is studied in 
higher dimensions. A study of the exact nature of the 
transition to backfiring and a comparison of the deter- 
ministic versus stochastic pathways to its existence will 
be undertaken in future work B. 

In summary, we proposed and studied a simple dis- 
crete model of calcium channel dynamics based on the 
assumption that calcium diffusion time is much smaller 
than the characteristic times of Ca 2+ binding/unbinding. 
This model demonstrates familiar properties of determin- 
istic reaction-diffusion systems in the limit N — > oo when 
fluctuations are small. For small N, we observed a tran- 
sition to a directed percolation regime, in agreement with 
the general DP conjecture p2| . For the full model includ- 
ing inhibition, we found at small N a novel persistent 
fluctuation driven state which emerges behind a front of 
outgoing activation; this occurs in a parameter regime 
where the corresponding deterministic system exhibits 
only single outgoing pulses. 



a 






' To whom correspondence should be addressed at 
falcke@mpipks-dresden.mpg.de 
[1] M. J. Berridge, M. D. Bootman, P. Lipp, Nature 395 645 

(1998) 
[2] J. Keizer, G.D.Smith, Biophys.Chem. 72 87 (1998) 
[3] J. Keizer, G.D.Smith, SPonce-Dawson, J.E.Pearson, 

1998, Biophys.J. 75 595 (1998) 
[4] M. Falcke, L. S. Tsimring, H. Levine, to be published. 
[5] G. W. DeYoung, J. Keizer, Proc. Natl. Acad. Set. USA 

89 9895 (1992) 
[6] J. Keizer, G. W. DeYoung, J. Theor. Biol. 166 431, 

(1994) 
[7] I. Bezprozvanny, J. Watras, B. E. Ehrlich, Nature, 351 

751 (1991) 
[8] X. P. Sun, N. Callamaras, J. S. Marchant, I. Parker, Jour- 
nal of Physiology 509.1, 67 (1998) 
[9] N. Callamaras, J. S. Marchant, X. P. Sun, I. Parker, Jour- 
nal of Physiology 509.1, 81 (1998) 
[10] E. T. Harris, Ann. Prob. 2 969 (1974) 
[11] S. R. Broadbent and J. M. Hammersley. Proc. Camb. 

Phil. Soc. 53, 629 (1957). 
[12] H. K. Janssen, Z. Phys. B 42, 151 (1981); P.Grassberger, 

Z.Phys. B 47, 365 (1982). 
[13] I. Mitkov, K. Kladko, and J. E. Pearson, Phys. Rev. Lett. 

81, 5453 (1999). 
[14] I. S. Aranson, B. Malomed, L. M. Pismen, L. S. Tsimring, 

submitted to Phys. Rev. Lett. . 
[15] M. G. Zimmerman, S.O. Firle. M. A. Natiello, M. Hilde- 
brand, M. Eiswirth, M. Bar, A. K. Bangia and I. G. 
Keverkides, Physica D 110, 92 (1997). 
[16] M. Bar, N. Gottschalk, M. Eiswirth and G. Ertl, J. 
Chern. Phys., 100, 1202 (1994). 



FIG. 5. Space-time evolution initiated by opening channels 
at a single cluster in the middle of the lattice of 300 sites 
for full activation/inhibition model with p + = 1, p^ — 0.04, 
p~ = 0.1, p^ = 0.12, h = 3, a = 0.7, and N s = 200 (a) and 
iV s = 20 (b), 500 iterations. 
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